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Abstract 

We study the interplay between traveling action potentials and spatial inhomogeneities in the 
FitzHugh-Nagumo model to investigate possible mechanisms for the occurrence of fibrillatory states 
in the atria of the heart. Different dynamical patterns such as ectopic foci, localized and meandering 
spiral waves are found depending on the characteristics of the inhomogeneities. Their appearance 
in dependence of the size and strength of the inhomogeneities is quantified by phase diagrams. 
Furthermore it is shown that regularly paced waves in a region R, that is connected by a small bridge 
connection to another region L with perturbing waves emanating from an additional pacemaker, 
can be strongly disturbed, so that a fibrillatory state emerges in region R after a transient time 
interval. This finding supports conjectures that fibrillatory states in the right atrium can be induced 
by self-excitatory pacemakers in the left atrium. 
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I. INTRODUCTION 



Atrial fibrillation (AF) is the most frequently appearing heart arrhythmia since it accounts 
'or one third of all hospitalizations caused by heart arrhythmia in the industrialized countries 
l|. During AF the electric conduction system of the heart is disturbed and an increased 
rate of activation by a factor of 3-12 compared to normal sinus rhythm occurs. Special 
spatio-temporal patterns of the electric potential like spiral waves^iother waves or ectopic 
foci are thought to be underlying generating mechanisms of AF [2-^. These patterns are 
often located near physiologically modified regions of the heart tissue in the left atrium js-Sj. 
The question hence arises, how these physiologically modified regions can be responsible for 
the generation of spiral waves or ectopic foci and how they influence the properties of these 
patterns. 

To tackle these questions, we study generating mechanism for AF on the basis of the 
FitzHugh-Nagumo model 9|, which is a simple model for action potential generation and 
propagation. By modeling physiologically modified regions using a spatial variation of the 
parameters characterizing cell properties like excitability or resting state stability, we calcu- 
late phase diagrams, which specify the type of spatio-temporal excitation pattern in depen- 
dence of the extent of the modified region and the strength of the modification. Thereupon 
we investigate how self-excitatory sources as spiral waves or ectopic foci with rather regular 
dynamics in one region can induce irregular, fibrillatory excitation patterns in some other 
region. Irregular, fibrillatory states are often observed in the right atrium , [lo| and it 
was conjectured that these are caused by the perturbation of regular waves generated by the 
sinus node by waves emanating from an additional pacemaker like a spiral wave or ectopic 
foci in the left atrium. 



II. MODEL 



The FitzHugh-Nagumo (FHN) equations 9| are a set of two coupled nonlinear ordinary 
differential equations, which describe excitable media via an inhibitor-activator mechanism. 
They were originally developed by searching for a simplified version of the Ho dgkin- Huxley 
equations for electric pulse propagation along nerves [llj. When combined with a spatial 
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diffusion term, the equations are 
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This set of partial differential equations serves as a prototype for a large variety of reaction- 
diffusion systems, which occur, for example, in chemical reactions as the Bhelousov- 



Zhabotinsky reaction 12l . Il3| or the catalysis of carbon monoxide 1^, ll5|], in population 



dynamics 
ics 



16| . in biology in connection with aggregation processes 17(| or plancton dynam 



191 ]. as well as in the spreading of forest fires 20 1. 



Here we will use Eqs. ([T]) in their original context as a model to investigate the spatio- 
temporal evolution of electric excitations in the heart. In this approach the variable u is 
roughly associated with the membrane potential and the variable v with the ion currents 
through the cell membrane. The resting state is given by the pair of values m = mq = 1.2 
and V = vq = —0.6. The diffusion coefficient D describes the coupling between the cells, and 
z is an applied external current (stimulus). The influence of the parameters a, b and c can 
be inferred by numerical solutions of Eqs. ([1]) without the diffusive term. The parameter 
values have to be limited to some range in order to generate excitability, and their detailed 
effect on the pulses is complicated due to mutual interdependencies originating from the 
nonlinearity in Eq. ([T]). Roughly speaking, a affects the length of the refractory period, b 
influences the stability of the resting state, and c controls the excitability and strength of 
the cells' response to a stimulus. To capture the propagation and form of a typical action 
potential, the following set of parameters can be used: D = Dq = 0.1, a = oq = 0.7, 
b = bo = 0.6, and c = cq = 5.5. These values will be associated with a "healthy tissue" in 
the following. Figure [1] shows the time development of u and v during an excitation after a 
stimulus with these parameters. Note that the variable —u mirrors the form of a pulse in 
an usual representation of an EGG recording. 

Ectopic foci and spiral waves are thought to be caused and influenced by physiologically 
modified regions of the tissue, which in the modeling correspond to spatial variations of the 
parameters. To simplify the analysis, we fix a = qq and D = Dq, and consider variations of 
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the parameters b and c according to 

b{x,y) = bo - Abexp{-^{x - xqY + {y - yo)y^b) , (2) 
c{x, y) = Co- Acexp(- {x - xoY + {y - yoY/ic) , (3) 

where the amphtudes A6, Ac characterize the strength, and the correlation lengths and 
characterize the spatial range of modification. 

The calculations are carried out on a two-dimensional simulation area of size 20 x 20, 
which represents an isolated section of atrial heart tissue, as it is used often in experiments 
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23| . The boundary conditions of the simulation area are of von Neumann type, i. e. 
du/dn = 0, where d/dn denotes the normal derivative. 

To solve the two nonlinear coupled partial differential equations ([T]) we use the finite 
element method (FEM) with a triangulation consisting of 4225 nodes and 8192 triangles, 
and a constant integration time step At = 0.01. A simulation time of 1 corresponds to a time 
of roughly 5 to 5.5 ms. The nonlinearity u^{x,t) in Eq. ([1]) is treated as an inhomogeneity, 
which means that for u{x,ti) the value u{x,ti_i) of the preceding time step is used. 

III. GENERATING MECHANISMS 
A. Ectopic activity 

Ectopic foci are regions in the atria, which generate activation waves emanating from 
self-excitatory hyperactive cells. In these cells the transmembrane potential raises without 
external stimulation until the threshold value is reached and an action potential results. 
In optical mapping studies and spatially resolved EGG recordings, ectopic foci are often 
localized in the regions of the pulmonary veins j?, sj. 

To model a tissue with physiologically modified properties that result in ectopic activity, 
we fix c = Co (Ac = 0) and vary the resting state stability b around the center of the 
simulation area (xq = yo = 10) according to Eq. ([2]). Initially the system is in the excitable 
resting state {u = Uo and v = Vo). Figure [2^ shows the resulting activation pattern for 
Ab = 0.4 and C,b = 0.8. The modified tissue is self-excitatory and acts as a pacemaker for 
activation waves, which propagate radially. In the time evolution shown in Fig. [3l u decreases 
until the threshold value Uth — 0.6 for activation is reached and an action potential with a 
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steep fall in u occurs. In response to this activation, the inhibitor variable v increases and 
pulls u back to a value even larger than its initial value (overshoot) before u returns to it, 
and the self-excitatory process starts anew. 

In order to systematically characterize the occurrence of ectopic activity, we calculate a 
phase diagram, where in dependence of C,b and Ab regions of ectopic activity can be distin- 
guished from that without self-excitatory behavior. The results in Fig. |2]d show that there 
exists a minimal A6min — 0.25, below which no ectopic activity occurs. The corresponding 
value (60 — A6niin) — 0.5 is the critical value of resting state stability in the FHN equations ([1]) 
in the absence of the diffusion term. For fixed Ab > Afemin the ectopic activity vanishes, 
when C,b falls below the dashed transition line in Fig. In this region the diffusive current 
from the modified tissue to the surrounding causes the initial decrease of u to become so 
slow, that the counter-regulation by v eventually hinders u to reach its activation threshold, 
cf. Fig. [3l Only small oscillation of u around a reduced resting state value can be seen in 
Fig. [31 which become weaker with growing time. 

The temporal-spatial pattern of the activation in the phase of ectopic activity is char- 
acterized by the frequency of the ectopic focus. We calculate this frequency as the inverse 
mean time interval between consecutive action potentials. As shown in Fig. HJ the frequency 
becomes larger with increasing Ab (at fixed ^f,) and C,b (at fixed Ab), and it tends to saturate 
for large With increasing C,b, the diffusive current of the inner cells of the modified tissue 
decreases and thus a larger frequency is obtained. In the saturation limit the frequency is 
nearly the same as in the absence of diffusion and thus is mainly determined by the refractory 
period. 



B. Spiral waves 



In this section we study the influence of physiologically modified regions, called "obsta- 
cles" henceforth, on spiral wave behavior. It was observed that spiral waves in the atria can 
be generated by a perturbation of the propagation of planar excitation waves by anatomic 
obstacles as, for example, the pulmonary veins, the venae cavae, the pectinate muscle bun 



dies or some localized region of modified tissue [6|, |21|, |24, |25| . These regions are considered 
as not fully excitable and are thus modeled as regions with a reduced parameter c according 
to Eq. 
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In an experiment by Ikeda and coworkers 2l| a nearly rectangular area of atrial tissue was 
placed on an electrode plaque in a tissue bath. Holes with different diameters were created 
and a reentrant wave was initiated by cross-field stimulation. The resulting behavior of the 
wavefront was, amongst others, classified according to whether the spiral is anchored by 
the obstacle, and by the relationship between hole size and cycle length of the reentry. It 
was observed that for large obstacle sizes (6, 8 and 10 mm), the reentrant wave attaches 
to the obstacle, leading to a linear increase of the cycle length with the hole diameter. For 
small obstacle diameters below about 4 mm, by contrast, meandering spirals with a tip 
getting variably closer to or further away from the hole were found. In this case the cycle 
length becomes independent of the hole diameter. Similar results were observed by Lim and 



coworkers 



24l |. They analyzed the behavior of spiral waves near holes with diameters ranging 



from 0.6 to 2.6 mm and obtained a higher attachment rate for larger obstacle diameters as 
well as a positive linear correlation of the reentry conduction velocity and wave length with 
the obstacle diameter in the case of attached spirals For smaller obstacle diameters the spiral 
waves were found to attach to and detach from the obstacle. 



The missing anchoring for small hole sizes was explained in (21| by invoking a "source-sink 
relationship". The "source" is the activation wavefront and provides a diffusive current to 
the surrounding tissue in the resting state, which constitutes the "sink". The sink becomes 
larger for smaller obstacles, where more cells become depolarized by the activation wavefront. 
If the source-to-sink ratio is decreased below a certain critical value, the wavefront detaches 
from the obstacle. 

To elucidate these experimental findings, we perform numerical calculations for a ge- 
ometry corresponding to the experiments with the following initial state and parameters 
settings: the modified region, is, as in the previous Sec. IIII At placed in the center of the 
simulation area at Xq = = 10. Initially a "planar" (linear) wave is generated by inducing 
a current z in the stripe 9.5 < a; < 10, < y < 10, and by setting the area < x < 9.5, 
< y < 10 into a refractory state with n = 1.6 and v = 0, while the rest of the simulation 
area is in the resting state {u = uq, v = vq). This initial state resembles the activation pat- 
tern directly after application of a cross-field or paired-pulse stimulation (two rectangular 
pulses). At the "upper part" of the initial planar wavefront (9.5 < x < 10, y = 10), diffusive 
currents flow "radially" in all forward directions {y > 10), while at the "right boundary" 
{x = 10, < y < 9.5) the diffusive currents can flow only in positive x direction (due to 



6 



the refractory state in the area < x < 9.5, < y < 10). This higher loss by diffusion 
leads to a smaller propagation speed of the initial wavefront at its upper boundary compared 
to its right boundary. As a consequence, the wavefront becomes curved, and a reentrant 
spiral wave develops for all reductions Ac in excitability and obstacle sizes in accordance 
with the experimental observations. Figure shows activation patterns for = 2 and a) 
Ac = 4.5, and b) Ac = 1.5. The stronger reduction of excitability in Fig. |5^ leads to an 
anchoring of the spiral wave, while in Fig. |5]d the spiral is meandering. 

To analyze the parameter regimes of the occurrence of anchored or meandering spiral 
waves, we perform a frequency analysis for different values of Ac and ^c- Therefore, we 
determine the peak positions in the time series of m at 8 positions far away from the center 
of the spiral and calculate the peak-to-peak intervals. The frequency of one point is one over 
the mean of the peak-to-peak intervals and the mean cycle length is one over the average of 
all these local frequencies. The results in Fig. [6] show that, as in the experiments, attached 
spiral waves occur for large Ac ^ 3 and for sufficiently large > where ^* decreases 
with increasing Ac. For these anchored spirals, the frequency is proportional to / = ri/27iC,c, 
where t] 0.82 is the conduction velocity in the FHN model. Accordingly, 1/f increases 
linearly with C,c for C,c > ■Cc Fig.[6l For small Ac ^ 3, only meandering spirals are observed. 
The transition from large to small Ac reflects the transition from anatomical to functional 



reentry 26|], as it has been reported in medical studies 2J]. The fact that for small always 
meandering spirals occur, can be interpreted by the small source-sink ratio 2l|. The same 
mechanism can also lead to meandering spirals for large C,c if Ac becomes small. Note, that 
nevertheless the spiral wave is not anchored to the obstacle, its movement is still influenced 
by the obstacle. 

IV. INDUCED FIBRILLATORY STATES IN THE RIGHT ATRIUM 

In previous studies on interactions of paced waves with self-excitatory waves, the influence 
of the pacing on a spiral wave was studied 
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with the aim to suggest a possible therapy 



to suppress fibrillation or tachycardia. The pacing was applied to the region, where the spiral 
wave was located. It was found that the pacing leads to an annihilation of the reentrant 



activity or to a shift of the spiral core 



28|-|31| 



Here we investigate the perturbation of regular paced waves from a source representing 
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the sinus node by waves emanating from an additional pacemaker located in a distant region. 
Electrocardiogram recordings and their frequency analyses show that regular excitation pat- 
terns are often observed in the left atrium, where additional pacemakers like spiral waves or 
ectopic foci are located Js, sl, and that at the same time irregular, fibrillatory-like states in 
the right atrium occur |7|. This led to the conjecture that fibrillatory states can be induced 
in the right atrium by self-excitatory pacemakers in the left atrium. In this connection it is 
important to better understand how a fibrillatory state can occur, if regular paced waves, as 
generated by the sinus node, are disturbed by additional pacemaker waves. To this end, we 
consider the waves to be located in spatially separated regions that are connected by a small 
region. To be specific, we choose a simulation area of size 21 x 10, which is divided into 
three regions (see Fig. [7]). The rectangular area L with < a; < 10, < ?/ < 10 representing 
the left atrium, the rectangular area R with 11 <x<21,0<y< 10 representing the right 
atrium, and the small bridge B with 10<x<ll,4<y<6 representing the connection 
between the atria. The grid used in the finite element calculations consists of in total 8871 
nodes and 17350 triangles. 

We focus on situations where the pacemaker in the region L is located far outside the 
left part of the simulation area, so that the resulting wavefronts become "planar" (linear). 
In the simulation they are generated by application of a stimulating current z = —1 with 
duration tz = I = lOOAt and a period l//pcrt in the region x < 0.5 and < y < 10. The 
activation waves representing the pacemaker in region R are generated by the application 
of a current z = —1 with duration tz = 1 and period l//pace in the region 11 < x < 21 and 
y < 0.5. 

The irregularity of the resulting patterns in region R is quantified by calculat- 
ing the Shannon entropy of the distribution of local activation frequencies for every 
grid point in R. To this end we divide the frequency range into A'^b bins of size 
A =min{exp [0.626 + 0.41n(A^g - 1)]~^ , 0.01 } ^ bsj and calculate the probabilities 

^,^ ..(/,</</,^A) _ 

of finding frequency / in bin /, where A'g is the total number of grid points. The normalized 
entropy then is given by 

s = = ^ (5) 

Smn.x In A^h 
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For a single frequency {pi = Si^i^), s = 0, while for a chaotic activation pattern with a uniform 
distribution {pi = 1/N\y), s = 1. 

For small perturbation frequencies (/pert ^ 0.1) the influence of the activation wavefronts 
from the additional pacemaker onto the sinus node waves is almost negligible. Small defor- 
mations of the linear wavefronts are observed, but the measured frequencies are close to the 
pacing frequency, and the overall spatiotemporal pattern in R is regular. 

With increasing perturbation frequency the spatiotemporal pattern in region R becomes 
more irregular and a breakup of the regularly paced waves can occur. Figure [8] shows the 
time evolution of the excitation patterns for a perturbation frequency /pert = 0.105 and a 
pacing frequency /pace = 0.091. For small times, the waves are only slightly perturbed and 
the pattern remains regular, as can be seen from the four consecutive snapshots in Fig. [H^. 
At a later time, however, the perturbation by the waves from region L results in a breakup 
of the waves close to the bridge B in region R, as can be seen from the four consecutive 
snapshots in Fig. [8b. The onset of this breakup was found to occur at a time t ^ 160. 

In order to investigate how the spatial irregularity is reflected in the time evolution, we 
consider three different points Pi = (11.49, 5.34), P2 = (11.33, 6.34) and P3 = (11.94, 6.15) in 
region R. The time evolution of u for these three points is shown for two different perturba- 
tion frequencies /p^t = 0.1 and /p^t = 0.105 in Fig. [9] (Pi: black sohd curve, P2: red dashed 
curve, P3: blue dash-dotted curve). For the lower frequency /p2t the evolution at all three 
points is regular, see Fig. [9^. For the higher frequency /pert? by contrast, the break ups of the 
waves seen in Fig. |8b yield unsuccessful activations during refractory periods, as can be seen, 
for example, at time t ^ 175 in point P2 (red dashed curve) and at time t = 179 in point 
Pi (black solid curve). These unsuccessful activations are caused by a rapid pacing of the 
region by the curled wave. Another feature is that the shape of the action potential varies. 
This can be seen, for example, at point Pi (black solid line) when comparing the pulses at 
t ^ 160 and t ~ 194. It is important to note that these irregularities are hardly observed at 
point P3 (blue dash-dotted curve), showing that they exhibit a spatial heterogeneity. 

The total irregularity in region R quantified by the normalized Shannon entropy s of the 
local frequency distribution is shown in Fig. [10] as a function of the frequency /pert of the 
perturbing waves from region L. For small perturbation frequencies /pert ^0.1 the entropy 
s equals the unperturbed case, while for /pert ^ 0.1, s sharply increases until it reaches a 
maximum at /pert — 0.1075. For higher /pert a return to more regular activation pattern is 
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found, indicating that the disturbance is most pronounced if /pert is close to /pace- 

To conclude, the disturbance of the wavefronts in region R by waves emanating from 
an additional source in region L and propagating through the bridge region B can lead to 
irregular, fibrillatory-like activation patterns in region R. On the other hand, the waves in 
region L are almost unaffected by the waves in R with primary wavefront orthogonal to the 
cross section of the bridge. The irregularities in region R are most pronounced at a certain 
perturbation frequency /pert- How this value is influenced by the geometry of the bridge and 
the parameters characterizing the cell properties remains further investigation. 



V. SUMMARY 



The influence of physiologically modified regions on the generation and properties of 
spatio-temporal activation patterns was investigated on the basis of the FitzHugh-Nagumo 
equations with von Neumann boundary conditions, in particular the occurrence of ectopic 
foci and spiral waves under spatial inhomogeneities of the parameters characterizing the 
cell properties. It was shown that the reduction Ab of the resting state stability in circular 
regions of the tissue can lead to ectopic activity. A minimal size of hyperactive tissue is 
necessary for ectopic activity to occur, as well as a minimal strength of the reduction of 
resting state stability with respect to the "healthy" reference value. With increasing size 
of the hyperactive tissue, the frequency of the ectopic focus first increases and eventually 
saturates. The saturation frequency depends on the strength of the modification Ab. 

For spiral wave patterns it was found that an anchoring of the wave to the obstacle can 
occur. To uncover this mechanism, an obstacle was modeled as a patch of modified tissue 
with reduced excitability by a reduction of the parameter c in the FHN equations. The 
obstacle was placed in the middle of a two-dimensional square simulation area and a planar 
excitation wave was generated aside of the obstacle in front of a refractory region, which 
represents an activation pattern observed after cross-field stimulation in experiments. As in 
the experiments, reentrant waves are observed. These exhibit either functional or anatomical 
reentry in dependence of the obstacle size and reduction strength Ac of excitability. An 
analysis of the spiral wave frequency in dependence of the obstacle size yields results in 
accordance with the experimental observations. 

Finally we studied the question, if and how fibrillatory-like states can arise in the right 
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atrium due to the presence of self-excitatory spiral waves or ectopic foci in the left atrium. To 
this end the simulation area was separated into two rectangular regions L and R connected 
by a small bridge B. Planar excitation waves were generated with different frequencies in 
the left region L to model a pacemaker far outside the left part of the simulation area. 
Planar excitation waves resembling stimulation by the sinus node were generated by periodic 
application of a stimulating current at one boundary in the right part R of the simulation 
area. For small perturbation frequencies in L, the disturbance of the waves in R turned out to 
be small. For higher perturbation frequencies , the waves in R become significantly disturbed 
and the spatio-temporal activation pattern eventually becomes irregular. The time evolution 
of the activation variable u, representing the electric potential in the FHN equations, shows 
features in close resemblance to the ones found in intra-atrial electrocardiogram recordings 
during fibrillation in the right atrium. The spatial variation of the excitation frequency was 
quantified in terms of an entropy, which showed, for a given pacing frequency, a maximum as 
a function of the perturbation frequency. Further investigations will focus on the influence 
of the geometry of the bridge and the wavefronts as well as analyse the behavior for different 
pacing frequencies. The reliability of the measure of irregularity s should be analysed and 
if necessary other methods to characterise the system behavior should be searched. 
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FIG. 1: Time evolution of a) u and b) v calculated with the FHN equations and the parameters 
representing a healthy tissue. 
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FIG. 2: (color online) a) Activation pattern of the ectopic focus at Ab = 0.4 and = 0.8. The 
lines represent isolines for u = —0.8 and four times t = 5 (black solid line), 7 (red dotted line), 9 
(blue dashed line), and 13 (green dash-dotted line). Initially {t = 0) the system is in the resting 
state {u = uo, V = vo)- b) Phase diagram of ectopic activity for modifications according to Eq. ([2]). 
Red diamonds and black circles refer to the occurrence and absence of ectopic activity, respectively. 
The dotted line is drawn as a guide to the eye and marks the transition between the regions of 
ectopic activity and absence of self-excitatory behavior. 
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FIG. 3: (color online) Time evolution of u and v for Ab = 0.4, and two values = 0.8 (black solid 
line) and ^j, = 0.5 (red dashed line). For = 0.8, u decreases below the threshold value Uth — 0.6 
and an action potential is generated. For ^f, = 0.5 the decrease of u is slower and less steep due to 
the stronger diffusive current compared to = 0.8. As a consequence the response of u is more 
susceptible to the initial decrease of v, which is not distinguishable for = 0.5 and = 0.8. Thus, 
u does not reach the threshold value and relaxes with small oscillations to a value u ~ 0.9. 
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FIG. 4: (color online) Frequency of ectopic activity in dependence of the size of the modified 
tissue for various modification strength Ab of the resting state stability. A frequency of 0.053 in 
the simulation corresponds to a frequency of roughly 10 Hz. 
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FIG. 5: a) Isolines for u = —0.8 at four different times t = 82 (black solid line), 87 (red dotted 
line), 92 (blue dashed line) and t = 97 (dash-dotted line) for an obstacle with = 2 and Ac = 4.5 
(marked by the black circle) . The spiral wave is pinned and rotates around the obstacle (anatomical 
reentry), b) Isolines for u = —0.8 at four different times t = 82 (black solid line), 84.5 (red dotted 
line), 87 (blue dashed line) and t = 89.5 (dash-dotted line) for an obstacle with (,c = 2 and Ac = 1.5. 
The spiral wave rotates around a moving center not corresponding to but influenced by the obstacle 
and is not attached to the obstacle (functional reentry). 
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FIG. 6: (color online) Mean cycle length 1/f as a function of obstacle size for four different 
strength of the modification of excitability Ac. A cycle length of 15 corresponds to a frequency of 
roughly 12 Hz. For the two smaller values Ac = 1.5 and Ac = 2.5 the cycle length is independent 
of (functional reentry). For the two larger values Ac = 3.5 and Ac = 4.5 a transition from 
functional to anatomical reentry occurs, when exceeds a threshold value ^* that increases with 
decreasing Ac. The dashed lines are fits to the data for mean cycle lengths independent of 
(functional reentry) and mean cycle lengths proportional to .^c (anatomical reentry). 
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FIG. 7: Illustration of the simulation area. The red dotted line describes a regular paced wavefront 
of the sinus node. The blue dashed line is a wavefront of the perturbing pacemaker. Pi, P2 and 
P3 mark the "observation points", for which the time evolution of u is shown in Figure [H 
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FIG. 8: Time evolution of excitation pattern observed by perturbation of a regular pacemaker with 
frequency /pace = 0.091 in the right part R of the simulation area by an additional pacemaker with 
frequency /pert = 0.105 in the left part L through a bridge region B. The blue solid lines are isolines 
for u = —0.8. The shaded regions represent the boundaries of the bridge between the two parts L 
and R. a) For small times the excitation in R is regular, b) At later times breakups of waves in 
region R occur close to the bridge region B, resulting in irregular excitation patterns. 

a) 2 = 




200 



FIG. 9: (color online) Time evolution of the activation variable u at three different points of the 
simulation area: Pi = (11.49,5.34) (black solid line), P2 = (11.33,6.34) (red dashed line) and 



P3 = (11.94,6.15) (blue dash-dotted line) for a pacing frequency /j 



paec 



0.091 and two different 



perturbation frequencies a) /pert = 0.1, and b) /pert = 0.105. We have checked that this behaviour 
remains qualitatively the same even at a three times longer simulation time which suggests that 
this time evolution corresponds to a stationary state. 
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FIG. 10: (color online) Normalized entropy of the local frequency distribution in region R as a 
function of the frequency /pert of the perturbing waves in region L. The pacing frequency in region 
R is /pace = 0.091. The dotted (black) line through the data points is a guide to the eye and 
the dashed line marks the value of s for the regularly paced system without perturbation by an 
additional pacemaker. 
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